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TESTS  OF  PARAMETERIZED  LANGMUIR-CIRCULATION  MIXING  IN 
THE  OCEAN’S  SURFACE  MIXED  LAYER 


1.  INTRODUCTION 

Martin  (1985,  1986)  reported  that  Mellor-Yamada-  (MY-)  type  turbulence  or  mixed-layer  mod¬ 
els  (MLMs)  (Mellor  and  Yamada  1974,  1982)  tend  to  underpredict  the  depth  of  the  surface  mixed 
layer  (SML)  in  tests  conducted  with  data  from  Ocean  Weathership  Stations  (OWSs)  November 
and  Papa.  Simulation  of  the  observed  response  of  upper-ocean  mixing  to  hurricanes,  e.g.,  Martin 
(1982),  also  suggested  that  these  types  of  MLMs  sometimes  do  not  mix  quickly  enough.  Observa¬ 
tions  of  velocity  shear  in  the  SML  have  suggested  that  the  velocity  shear  in  the  SML  predicted  by 
MY-type  models  is  too  large,  which  also  suggests  that  the  rate  of  mixing  is  not  strong  enough. 

There  have  been  a  number  of  efforts  to  enhance  the  MY-type  MLMs  to  account  for  the  fact 
that  they  frequently  do  not  mix  the  ocean’s  SML  strongly  or  deeply  enough.  Large  et  al.  (1994) 
proposed  a  mixing  enhancement  to  account  for  possible  intermittent  mixing  in  the  stratified  region 
near  the  base  of  the  SML  by,  e.g.,  enhancement  of  the  local  shear  by  propagating  or  breaking 
internal  waves.  This  enhancement  provides  for  a  moderate  rate  of  vertical  mixing  (with  a  maximum 
rate  of  50  cm2/s)  up  to  a  Richardson  number  of  0.7,  i.e.,  the  additional  mixing  provided  by  the 
enhancement  shuts  off  at  a  Richardson  number  of  0.7,  rather  than  a  value  in  the  range  of  0.20  to  0.25 
that  is  typical  of  most  MY-type  MLMs.  This  mixing  enhancement  was  implemented  and  tested  in 
their  K-Profile  Parameterization  (KPP)  turbulence  model  by  Large  et  al.  (1994)  and  was  found  to 
significantly  improve  the  predicted  mixed-layer  depth  (MLD)  relative  to  observations.  The  Large 
et  al.  (1994)  mixing  enhancement  was  also  implemented  and  tested  in  the  Mellor-Yamada  Level 
2.5  (MYL2.5)  MLM  by  Kantha  and  Clayson  (1994),  and  they  reported  a  significant  improvement 
in  the  predicted  MLD  for  several  data  sets  when  using  this  enhancement. 

Another  enhancement  to  increase  the  mixing  for  MY-type  MLMs  has  been  to  account  for  the 
injection  of  turbulent  kinetic  energy  (TKE)  into  the  SML  by  breaking  surface  waves.  This  can  be 
parameterized  in  the  MYL2.5  and  similar  turbulence  models  by  implementing  a  flux  of  TKE  at  the 
surface  from  the  waves  and  by  setting  the  turbulent  length  scale  at  the  surface  zs  to  account  for 
the  surface  roughness  due  to  the  surface  waves. 

Craig  and  Banner  (1994)  parameterized  the  surface  flux  of  TKE  as  cm*,  where  a  =  100  and 
u*  is  the  water-side  surface  friction  velocity.  For  their  experiments  with  u*  =  0.011  m/s,  they  used 
a  value  of  0.1  m  for  the  surface  roughness.  However,  they  note  that  this  choice  was  somewhat 
arbitrary  and  that  there  is  considerable  uncertainty  in  defining  the  surface  roughness  for  upper- 
ocean  turbulence,  i.e.,  their  discussion  of  estimates  from  the  literature  for  light  to  moderate  winds 
includes  values  from  10-4  to  more  than  1  m.  They  note  that  the  most  popular  approach  to 
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computing  zs  is  probably  the  Charnock  (1955)  formula  zs  =  f3u*/g,  where  / 3  is  a  scaling  factor  and 
g  is  the  acceleration  of  gravity. 

The  initial  finding  of  using  this  approach  to  increasing  the  mixing  in  MY-type  MLMs  is  that 
the  additional  TKE  added  near  the  surface  increases  the  mixing  rates  near  the  surface,  but  that  this 
additional  TKE  tends  to  dissipate  locally,  i.e. ,  near  the  surface,  and  does  not  generally  contribute 
significantly  to  increasing  the  vertical  mixing  deeper  in  the  water  column  or  increasing  the  depth 
of  the  SML  (Craig  and  Banner  1994).  However,  if  the  surface  roughness  is  made  sufficiently  large 
(i.e.,  larger  than  most  of  the  values  that  were  initially  tested),  the  mixing  rates  throughout  the 
SML  are  significantly  increased  and  the  depth  of  the  SML  is  increased  as  well. 

This  issue  was  reviewed  by  Mellor  and  Blumberg  (2004).  Based  on  predicting  a  sufficiently 
deep  SML  at  OWS  Papa  with  the  MYL2.5  MLM,  they  recommended  values  of  a  =  100  for  the 
TKE  flux  (i.e.,  the  same  as  the  value  suggested  by  Craig  and  Banner,  1994)  and  (5  =  200000  for 
the  calculation  of  the  surface  roughness  using  the  Charnok  (1955)  formula.  Note  that  this  value  of 
/ 3  gives  a  surface  roughness  of  about  2  m  in  winds  of  8  m/s,  which  is  of  the  order  of  the  surface 
wave  height.  Whether  such  a  large  surface  roughness  is  appropriate  for  computing  turbulence  in 
the  SML  is  not  clear.  However,  these  large  values  of  [3  do  provide  a  simple  way  to  increase  the 
depth  of  the  SML  predicted  by  MY-type  turbulence  models,  and  the  predicted  deepening  of  the 
SML  can  be  adjusted  by  adjusting  the  value  of  /3.  Note  that  for  the  simulations  conducted  in  this 
report,  the  value  of  a  was  set  to  100  and  the  value  of  j3  was  set  to  40000.  This  value  of  j3  is  larger 
than  the  smallest  values  that  have  been  suggested,  but  is  smaller  than  the  largest  and  is  smaller 
than  is  needed  to  significantly  affect  the  MLD. 

It  has  long  been  thought  that  Langmuir  circulations  (LCs)  must  affect  upper-ocean  mixing. 
However,  until  fairly  recently,  there  have  not  been  quantitative  estimates  as  to  what  the  effects  of 
LCs  on  upper-ocean  mixing  are,  and  so  there  have  been  few,  if  any,  attempts  to  parameterize  the 
effect  of  LCs  on  vertical  mixing  in  the  SML. 

However,  improved  understanding  of  LCs  and  increases  in  computing  power  have  led  to  large- 
eddy  simulations  (LES)  of  LCs,  which  have  provided  quantitative  estimates  of  the  effect  of  LCs  on 
mixing  in  the  ocean’s  SML.  McWilliams  et  al.  (1997,  hereafter  referred  to  as  MW97)  conducted 
LES  simulations  of  LCs  for  a  simple,  wind-mixing  case  with  an  initial  MLD  of  about  33  m  and  a 
moderate  wind  of  about  5  m/s.  They  found  that  the  inclusion  of  LCs  in  the  LES  increased  the 
rate  of  mixing  within  the  SML  by  about  a  factor  of  3  and  slightly  increased  the  depth  of  mixing. 

Motivated  by  this  result  from  MW97,  Kantha  and  Clayson  (2004,  hereafter  referred  to  as 
KC04)  investigated  how  to  implement  this  increased  rate  of  mixing  from  the  LCs  in  the  widely- 
used  MYL2.5  turbulence  model.  They  added  additional  shear-production  terms  to  the  MYL2.5 
model,  consisting  of  the  product  of  the  vertical  Reynolds  stress  multiplied  by  the  vertical  shear  of 
the  Stokes  drift  current  (SDC)  from  the  waves.  The  additional  shear-production  terms  were  added 
to  both  the  TKE  equation  and  the  turbulence  length-scale  (TLS)  equation  of  the  MYL2.5  model. 

The  effect  of  these  changes  on  the  MYL2.5  turbulence  model  is  generally  to  increase  the  TKE 
and  the  TLS  and  the  rate  of  mixing  within  the  SML.  For  the  conditions  of  the  simple,  light-wind 
simulation  conducted  by  MW97,  KC04  found  that  the  maximum  rate  of  mixing  in  the  SML  was 
significantly  increased  and  the  mixing  rates  were  fairly  consistent  with  the  LES  results  of  MW97. 
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Some  properties  of  the  KC04  modification  to  the  MYL2.5  turbulence  model  are  that  (1)  since 
it  depends  on  the  SDC  from  the  waves,  the  additional  mixing  is  generally  confined  to  the  region 
of  influence  of  the  surface  waves  in  the  SML,  i.e. ,  the  additional  mixing  does  not  affect  interior  or 
bottom  mixed  layers  unless  the  water  is  fairly  shallow,  and  (2)  since  the  additional  shear- production 
terms  depend  on  the  product  of  the  vertical  momentum  flux  and  the  vertical  shear  of  the  SDC,  and 
the  vertical  momentum  flux  is  parameterized  as  the  product  of  the  vertical  mixing  coefficient  and 
the  vertical  shear  of  the  ocean  model’s  Eulerian  velocity,  if  there  are  waves  but  the  winds  are  light, 
so  that  the  vertical  shear  of  the  model  Eulerian  velocity  is  small,  the  additional  shear-production 
terms  will  not  generate  much  mixing,  i.e.,  in  the  presence  of  waves  or  swell  but  weak  winds,  the 
additional  shear-production  terms  will  generally  not  cause  much  mixing.  Both  of  these  properties 
may  be  desirable  for  agreement  with  observations  of  upper-ocean  mixing. 

The  modifications  to  the  MYL2.5  turbulence  model  prescribed  by  KC04  were  implemented 
in  the  version  of  the  MYL2.5  turbulence  model  used  in  the  Navy  Coastal  Ocean  Model  (NCOM) 
(Martin  2000;  Morey  et  al.  2003;  Barron  et  al.  2004)  and  some  tests  were  conducted  to  see  how  the 
prescribed  changes  perform.  The  following  sections  contain  a  description  of  the  implementation  of 
the  KC04  LC  mixing  enhancement  in  the  MYL2.5  turbulence  model  in  NCOM  (Section  2),  a  test 
of  the  KC04  LC  mixing  in  NCOM  for  the  simple,  light-wind,  test  case  used  by  MW97  and  KC04 
(Section  3),  a  test  of  the  KC04  LC  mixing  in  NCOM  at  OWS  Papa  (Section  4),  a  test  of  the  KC04 
LC  mixing  in  NCOM  for  a  simulation  of  Hurricane  Ivan  in  the  Gulf  of  Mexico  (Section  5),  and  a 
summary  (Section  6). 

Note  that  the  goal  of  the  tests  is  to  see  how  much  difference  the  KC04  LC  mixing  parameter¬ 
ization  makes  relative  to  using  the  basic  form  of  the  MYL2.5  turbulence  model.  In  this  way  we 
generally  follow  the  testing  procedure  used  by  KC04,  who  chose  not  make  use  of  the  internal-wave 
mixing  enhancement  of  Large  et  al.  1994,  which  they  had  previously  implemented  in  the  MYL2.5 
MLM  to  increase  the  depth  of  mixing  (Kantha  and  Clayson  1994).  Hence,  underprediction  of  the 
observed  MLD  in  the  tests  conducted  for  this  report,  e.g.,  at  OWS  Papa,  is  not  to  be  considered 
a  significant  issue,  since  such  shortcomings  can  be  remedied  by  incorporating  additional  mixing 
mechanisms,  e.g.,  the  internal- wave  mixing  enhancement  of  Large  et  al.  1994  or  the  large  surface 
roughness  suggested  by  Mellor  and  Blumberg  (2004)  or  some  other  mixing  mechanism. 

2.  IMPLEMENTATION  OF  LC  MIXING  IN  MYL2.5  MLM 


The  original  TKE  and  TLS  equations  as  implemented  in  NCOM  are 


dq2 

8t 


-v  •  (W)  +  Qq2  +  vh(AHvhq2)  +  ^ 


+  2  Km1 


+  2KHt<Lf 

Po  OZ 


dq2 1 

~df 


-V  •  (vq2l)  +  Qq2i  +  Vh{AHVh{q2l))  +  ^  (xqi 

+  EitKMi  ((£)’+  (|)2)  +E,eKHlf^-E,lW, 


(1) 


(2) 
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where  q  is  the  square  root  of  twice  the  TKE,  i  is  the  turbulent  length  scale,  t  is  the  time,  v  is 
the  vector  velocity,  Q  is  a  volume  flux  source  term,  Ajj  is  the  horizontal  mixing  coefficient  for  the 
scalar  fields,  z  is  the  vertical  coordinate,  Kqi  is  the  vertical  diffusion  coefficient  for  q 2  and  q2i, 
which  is  taken  to  be  proportional  to  Km1  ( Kqt  =  0.41^^),  Km i  is  the  vertical  mixing  coefficient 
for  momentum,  u  and  v  are  the  horizontal  components  of  the  velocity,  Kji1  is  the  vertical  diffusion 
coefficient  for  the  scalar  fields,  g  is  the  acceleration  of  gravity,  p0  is  a  reference  density  for  the  water, 
and  dp/dz  is  the  vertical  gradient  of  the  water  density. 


W  is  referred  to  by  Mellor  and  Yamada  (1982)  as  a  “wall  proximity”  function,  which  is  used 
to  scale  l  near  the  surface  and  bottom.  This  function  is  defined  by 

W  =  1+Ei(Y)\  ,3) 

and  in  NCOM,  L  is  defined  by 

L  1  =  (C  —  z  +  zs)  1  +  (z  —  H  +  Zb)  1,  (4) 


where  n  =  0.4  is  Von  Karman’s  constant,  (  is  the  elevation  of  the  free  surface,  zs  is  the  surface 
roughness,  Zf,  is  the  bottom  roughness,  and  H  is  the  static  bottom  depth.  The  symbols  b\ .  E\ .  E2 , 
E%,  and  E4  are  constants  and  the  values  of  these  are  specified  in  Table  1. 

Table  1  —  Constants  for  MYL2.5  Turbulence  Equations 


parameter 

value 

bi 

16.6 

Ei 

1.8 

E2 

1.0 

E3 

1.8 

Ei 

1.33 

Eq 

7.2 

KC04  parameterized  the  effects  of  the  LC  mixing  reported  by  McWilliams  et  al.  (1997)  in  their 
LES  experiment  by  adding  additional  shear  production  terms  to  both  the  TKE  and  TLS  equations 
in  the  MYL2.5  turbulence  model.  The  TKE  and  TLS  equations  in  NCOM  with  the  additional 
shear  production  terms  used  by  KC04  are 
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where  us  and  vs  are  the  horizontal  components  of  the  SDC  from  the  waves.  An  additional  constant 
Ef,  is  used  to  scale  the  Stokes  shear  production  term  in  the  TLS  equation  and  has  a  value  of  7.2 
(Table  1).  Note  that  the  value  of  7.2  for  Eq  is  different  from  the  value  of  4.0  reported  in  KC04, 
which  was  incorrect  (Dr.  Lakshmi  Kantha,  personal  communication).  We  found  in  early  testing 
that  the  KC04  LC  mixing  parameterization  does  not  work  correctly  (i.e. ,  as  reported  by  KC04) 
with  Ef,  =  4.0. 

The  original  shear-production  terms  are  proportional  to  the  square  of  the  vertical  velocity  shear; 
hence,  the  original  shear-production  terms  are  always  positive  and  always  act  to  increase  the  TKE. 
The  shear-production  terms  added  by  KC04  to  parameterize  the  LC  mixing  contain  the  vertical 
shear  of  the  model  Eulerian  velocity  times  the  vertical  shear  of  the  SDC;  hence,  it  is  possible  for 
this  term  to  be  negative  if  the  two  velocity  shears  are  of  different  sign,  which  would  result  in  a  loss 
of  TKE  from  this  term,  rather  than  a  gain.  However,  in  the  tests  we  have  conducted  so  far,  this 
has  not  seemed  to  be  a  large  problem,  probably  because  the  vertical  shear  of  the  model  Eulerian 
current  and  the  SDC  usually  have  the  same  sign.  There  will  be  some  further  comments  about  this 
issue  in  the  test  results. 

The  coupling  of  NCOM  with  the  output  from  a  wave  model  requires  some  additional  changes 
besides  enhancing  the  vertical  mixing  within  the  MYL2.5  turbulence  model.  The  SDC  is  added  to 
NCOM’s  Coriolis  term  as  described  by  MW97  and  KC04  and  many  others.  This  term  is  referred  to 
as  the  Stokes-Coriolis  term  and  it  affects  the  downwind  transport  in  the  surface  wind-driven  layer 
as  discussed  in  MW97  and  KC04  and  in  this  report  in  the  next  section.  The  SDC  is  also  used  to 
advect  the  model  fields  in  NCOM  and,  to  be  consistent  with  this  advection,  is  added  to  NCOM’s 
continuity  equation.  Note  that  the  addition  of  the  SDC  to  NCOM’s  Coriolis  terms  affects  local, 
one-dimensional  simulations,  such  as  those  presented  in  this  report,  but  the  advection  by  the  SDC 
does  not. 

Additional  wave  effects  include  the  surface  wave-radiation  stress,  which  is  implemented  in 
NCOM  similar  to  the  surface  wind  stress;  but  since  it  is  generally  much  smaller  than  the  wind 
stress,  its  effect  is  usually  small.  Another  effect  of  the  waves  is  the  enhancement  of  bottom  drag  in 
shallow  water  due  to  the  surface-wave  orbital  motions  near  the  bottom  (Grant  and  Madsen  1979; 
Soulsby  1995).  Both  of  these  additional  wave-forcing  effects  are  implemented  in  NCOM,  but  these 
effects  are  not  the  focus  of  this  study  and  will  not  be  discussed  further  in  this  report. 

3.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  FOR  SIMPLE  WIND  MIXING  CASE 

MW97  used  a  simple,  wind-mixing  test  case  to  look  at  the  effects  of  LCs  on  vertical  mixing,  and 
KC04  ran  the  same  case  to  test  their  parameterization  of  LC  mixing  in  the  MYL2.5  MLM  and  to 
compare  their  results  with  those  of  MW97.  Hence,  we  ran  the  same  case  to  test  our  implementation 
of  KC04’s  parameterization  of  LC  mixing  in  the  version  of  the  MYL2.5  MLM  in  NCOM  and  to 
compare  with  the  results  of  MW97  and  KC04. 

The  initial  condition  for  this  test  case  consists  of  a  SML  with  a  depth  of  33  m  and  a  temperature 
of  13.5  °C  and  a  thermal  stratification  below  the  SML  of  0.01  °C/m.  The  salinity  is  set  to  a  uniform 
35  psu  both  within  and  below  the  SML.  The  initial  temperature  was  chosen  to  give  a  thermal 
expansion  coefficient  (1  / pdp/8T)  of  -0.0002  1/°C.  The  latitude  is  43.4  °N  to  give  a  value  of  the 
Coriolis  parameter  of  0.0001  1/s,  which  results  in  an  inertial  period  of  17.45  h. 
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The  forcing  for  the  test  case  consists  of  a  surface  wind  stress  of  0.037  Pa,  which  corresponds  to 
a  wind  speed  of  about  5  m/s.  MW97  used  a  small  cooling  surface  heat  flux  of  5  W/m2  to  speed  up 
the  convergence  of  their  LES  numerical  solution  and  reduce  its  run  time.  We  used  a  linear  ramp 
the  length  of  the  inertial  period  (17.45  h)  for  the  surface  wind  stress  to  reduce  inertial  oscillations 
and  speed  the  convergence  to  an  approximately  steady  solution. 

The  surface  wave  field  for  this  test  case  was  specified  by  MY97  as  a  monochromatic  wave  of 
amplitude  0.8  m  and  wavelength  60  m  propagating  in  the  downwind  direction.  This  gives  a  SDC 
with  a  surface  value  of  0.0679  m/s  and  an  e-folding  depth  (the  inverse  of  twice  the  wavenumber, 
see  Appendix  A)  of  4.78  m,  which  corresponds  to  a  light  swell.  This  is  the  SDC  used  by  KC04 
and  by  us  for  this  test. 

NCOM  was  run  for  this  case  in  what  we  refer  to  as  pseudo  one-dimensional  mode  by  running 
the  model  on  a  horizontal  domain  of  2  by  2  grid  points  (the  minimum  horizontal  grid  that  NCOM 
allows)  with  doubly-periodic  lateral  boundary  conditions.  The  vertical  grid  consisted  of  40  layers 
with  an  upper-layer  thickness  of  1  m  and  a  uniform  stretching  to  a  maximum  depth  of  200  m; 
hence,  each  layer  is  thicker  than  the  layer  above  it  by  about  7%. 

With  the  use  of  a  linear  ramp  for  the  surface  wind  stress,  the  numerical  solution  is  fairly  steady 
once  the  ramp  ends  after  17.45  h.  Figure  1  shows  profiles  of  TKE,  the  vertical  mixing  coefficient 
for  momentum,  and  the  downwind  and  crosswind  velocity  for  simulations  without  and  with  the 
parameterization  of  LC  mixing.  This  figure  can  be  compared  with  Fig.  2,  3,  and  4  in  MW97  and 
with  Fig.  1  in  KC04. 

The  TKE  in  Fig.  1  is  roughly  doubled  by  the  use  of  the  KC04  LC  mixing  parameterization, 
which  is  qualitatively  similar  to  the  result  in  KC04.  Note  that  the  increase  in  TKE  in  the  upper 
meter  in  Fig.  1  is  due  to  the  parameterization  of  the  surface  flux  of  TKE  from  the  surface  waves. 

The  maximum  value  of  the  vertical  mixing  coefficient  for  momentum  Km  in  Fig.  1  is  also 
roughly  doubled,  which  is  qualitatively  similar  to  the  result  in  KC04,  i.e. ,  the  maximum  value  of 
Km  increases  from  about  200  to  400  cm2/s,  compared  with  an  increase  from  about  215  to  430  cm2/s 
for  KC04.  For  MY97  (their  Fig.  3b),  the  maximum  value  of  Km  increases  by  about  a  factor  of  3 
from  about  130  to  380  cm2/s.  However,  their  maximum  value  with  the  LC  mixing  (380  cm2/s)  is 
similar  to  that  obtained  in  our  simulation  (400  cm2/s)  and  by  KC04  (430  cm2/s).  The  shape  of  the 
vertical  mixing  coefficient  profiles  in  Fig.  1  is  more  similar  to  those  in  MW97  than  in  KC04  in  that 
the  profiles  in  Fig.  1  and  in  MW97  have  a  roughly  parabolic  shape,  whereas  the  profiles  in  KC04 
reduce  fairly  abruptly  from  a  maximum  value  to  a  value  close  to  zero  within  a  short  distance  near 
the  base  of  the  SML. 

The  downwind  velocity  u  in  Fig.  1  shows  reduced  shear  with  the  stronger  LC  mixing  and,  as 
a  result,  the  value  of  u  at  the  surface  is  noticeably  reduced,  which  is  also  shown  by  KC04  (their 
Fig.  1)  and  MW97  (their  Fig.  2).  Note  that  the  velocity  profiles  in  Fig.  1,  in  Fig.  2  in  MW97,  and 
in  Fig.  1  in  KC04  are  for  just  the  models’  Eulerian  current  and  do  not  include  the  SDC. 

With  the  inclusion  of  the  SDC  in  the  Coriolis  term  in  the  momentum  equations,  the  downwind 
transport  of  the  Eulerian  current  is  not  zero  but  balances  the  transport  of  the  SDC,  so  that  the 
total  downwind  transport  of  the  Eulerian  current  plus  the  SDC  is  zero.  Hence,  since  the  transport 
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V 


V  (cm/s) 


Fig.  1  —  Profiles  of  TKE,  vertical  mixing  rate  for  momentum,  and  downwind  and  crosswind  velocity  for  simulations 
without  (solid  line)  and  with  (dashed  line)  the  KC04  LC  mixing  parameterization  for  the  simple,  light-wind,  test 
case. 
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of  the  SDC  is  downwind  for  this  particular  test  problem,  the  transport  of  the  downwind  component 
of  the  Eulerian  current  in  the  models  is  upwind  as  shown  in  Fig.  1  and  by  KC04  and  MW97. 

A  difference  of  MW97  from  the  simulation  here  and  in  KC04  for  the  case  without  the  LC 
mixing  is  that  MW97  drop  the  SDC  from  the  simulation  completely,  whereas  here  and  in  KC04  the 
SDC  is  retained  in  the  calculation  of  the  Coriolis  term  and  is  only  turned  off  for  the  calculation  of 
the  vertical  mixing  in  the  MYL2.5  turbulence  model.  Hence,  in  Fig.  1  and  in  KC04’s  Fig.  1,  the 
transport  of  u  for  the  case  without  LC  mixing  is  still  upwind  to  balance  the  downwind  transport 
of  the  SDC,  whereas  in  MW97,  the  transport  of  u  for  the  case  without  LC  mixing  is  zero  (their 
Fig.  3). 

The  cross-wind  velocity  v  in  Fig.  1  also  shows  reduced  shear  with  the  LC  mixing,  similar  to 
the  cross-wind  velocity  in  MY97  (their  Fig.  2),  and  the  surface  value  of  the  cross-wind  velocity  is, 
as  a  result,  reduced.  The  transport  of  the  cross-wind  velocity  must  equal  the  Ekman  transport  for 
all  of  these  simulations  once  they  have  reached  near-equilibrium. 

The  profiles  in  Fig.  1  indicate  slightly  deeper  mixing  for  the  case  with  the  LC  mixing,  i.e.,  for 
the  case  without  the  LC  mixing  the  MLD  increases  from  the  initial  value  of  33  to  34  m,  and  for  the 
case  with  the  LC  mixing  the  MLD  increases  from  33  to  37  m.  A  deeper  MLD  with  the  LC  mixing 
is  also  indicated  in  the  results  of  MW97  and  KC04. 

The  velocity  profiles  in  KC04  (their  Fig.  1)  do  not  appear  completely  consistent,  since  the 
transports  for  the  cases  with  and  without  LC  mixing  should  be  the  same  and  they  are  somewhat 
different.  This  may  be  because  the  velocity  profiles  shown  by  KC04  have  not  reached  equilibrium 
and  there  is  some  residual  inertial  motion  present  in  the  velocity  profiles  that  they  show. 

In  summary,  the  results  obtained  for  this  simple,  wind-mixing  test  case  are  generally  similar  to 
the  results  reported  by  KC04  and  MW97.  The  values  of  the  maximum  mixing  rate  for  momentum 
are  roughly  doubled  by  including  the  LC  mixing  in  the  results  here  and  in  KC04.  This  is  less  than 
the  factor  of  three  increase  obtained  by  MW97,  but  the  maximum  mixing  rate  for  momentum  with 
the  LC  mixing  obtained  here  of  400  cm2/s  is  similar  to  that  obtained  by  KC04  and  MW97  with  the 
LC  mixing.  The  approximately  parabolic  shape  of  the  eddy  coefficient  mixing  profiles  obtained  here 
with  NCOM  are  roughly  similar  to  the  shape  of  the  eddy  coefficient  profiles  obtained  by  MW97, 
but  are  somewhat  different  from  the  shape  of  the  profiles  in  KC04. 

4.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  AT  OWS  PAPA 

OWS  Papa  was  located  in  the  northeast  Pacific  at  145  °W,  50  °N.  The  availability  of  long  time 
series  of  meteorological  and  ocean  subsurface  measurements  at  Papa,  and  the  fact  that  the  effects 
of  advection  on  the  heat  budget  tend  to  be  small  in  this  area,  have  long  made  Papa  a  popular 
location  to  test  and  evaluate  upper-ocean  MLMs  (Denman  and  Miyake  1973;  Mellor  and  Durbin 
1975;  Martin  1985;  Martin  1986;  Gaspar  1988;  Large  et  al.  1994;  Kantha  and  Clayson  1994;  Large 
1996). 

Simulations  were  conducted  at  OWS  Papa  for  the  year  1961  using  NCOM  run  in  pseudo  one¬ 
dimensional  mode  and  using  the  MYL2.5  turbulence  model  both  with  and  without  the  LC  mixing 
parameterization  of  KC04  (i.e.,  similar  to  the  simulations  discussed  in  the  previous  section).  Hourly 
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surface  wind  stresses  and  heat  fluxes  were  computed  from  the  3-hourly  surface  marine  observations 
at  Papa  using  fairly  standard  formulas  as  discussed  in  Martin  (1985,  1986). 

Note  that  all  the  heat  fluxes  were  pre-computed  using  the  surface  marine  observations  from 
Papa,  i.e. ,  the  NCOM-predicted  SST  was  not  used  in  the  calculation  of  any  of  the  heat  fluxes.  The 
mean  computed  heat  fluxes  for  the  year  were  105  W/m2  for  the  solar  radiation,  -88  W/m2  for  the 
surface  heat  loss  due  to  the  net  longwave  and  latent  and  sensible  heat  fluxes,  and  18  W/m2  for  the 
net  heat  flux.  These  pre-computed  heat  fluxes  are  fairly  consistent  with  the  seasonal  changes  in 
the  observed  upper-ocean  heat  content  at  Papa  during  the  spring  and  summer  of  1961  (Martin  et 
al.  2012). 

Solar  extinction  was  taken  to  be  Jerlov  Type  II  (Jerlov  1968)  and  the  ambient/background 
viscosity  and  diffusivity  below  the  SML  were  taken  to  be  0.02  cm2/s.  A  linear  damping  term  in  the 
momentum  equations  with  a  damping  time  scale  of  10  d  was  used  to  provide  some  damping  of  the 
inertial  oscillations,  since  the  normal  horizontal  and  vertical  dissipation  of  these  motions  cannot  be 
accounted  for  in  a  one-dimensional  simulation  such  as  is  being  conducted  here. 

The  initial  condition  for  the  simulations  at  Papa  was  a  temperature  profile  from  the  bathyther¬ 
mograph  (BT)  observations  at  Papa  for  1  January  1961,  and  a  climatological  salinity  profile  for  the 
location  of  Papa  for  the  beginning  of  January  from  the  Levitus  (1982)  climatology.  The  salinity 
profiles  at  Papa  feature  a  strong  halocline  between  about  100  and  200  m  that  helps  limit  the  depth 
of  mixing  in  the  winter. 

Including  the  effect  of  the  LC-mixing  parameterization  of  KC04  on  the  simulation  of  the  SML 
at  Papa  requires  an  estimate  of  the  SDC.  Since  we  don’t  have  wave  observations  for  Papa  for  1961, 
the  SDC  profile  was  estimated  from  the  observed  surface  winds  at  Papa.  There  are  a  number  of 
ways  this  can  be  done,  and  three  methods  are  compared  here:  (1)  compute  the  SDC  from  the  wind 
stress  using  the  method  used  by  KC04  for  their  simulations  at  OWS  Papa,  (2)  compute  the  SDC 
from  the  wind  speed  using  wave  energy  spectra  time-averaged  for  different  wind  speeds  using  wind 
and  wave  data  from  a  NOAA  Buoy  located  in  the  northeast  Pacific  about  1100  km  east-southeast  of 
Papa,  and  (3)  compute  the  SDC  from  the  waves  predicted  by  a  wave  model  forced  by  the  observed 
winds  at  Papa  for  1961. 

4.1  Papa  Simulations  with  SDC  Computed  Using  Method  of  KC04 

The  method  used  to  estimate  the  SDC  from  the  surface  wind  stress  by  KC04  for  their  simu¬ 
lations  at  OWS  Papa  to  test  their  LC  mixing  parameterization  is  to  define  the  surface  magnitude 
of  the  SDC  Vs  as  a  linear  function  of  the  friction  velocity  u*,  i.e.,  14,(0)  =  11. 8w*,  and  the  depth- 
dependence  of  the  SDC  is  defined  by  the  wave  number  k  that  is  characteristic  of  the  peak  of  the 
wave-energy  spectrum,  which  is  parameterized  as  k  =  4.05  x  10 ~6g/ul  (see  Appendix  B).  The 
direction  of  the  SDC  is  taken  to  be  downwind.  This  will  be  referred  to  as  the  KC04  SDC.  Using 
this  method,  hourly  profiles  of  the  SDC  were  computed  from  the  surface  wind  stress  at  Papa  for 
1961  and  stored  in  a  file  for  input  to  NCOM. 

Figure  2  shows  results  from  a  simulation  of  the  SML  at  Papa  for  1961  using  just  the  MYL2.5 
MLM  with  no  LC  mixing  enhancement.  However,  note  that  the  SDC  is  used  in  the  Coriolis  term 
of  the  ocean  model  as  discussed  in  Section  3.  In  the  figure,  the  results  of  the  simulation  (red)  are 
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Fig.  2  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  KC04  SDC,  but  without 
the  KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 


compared  with  the  observations  at  Papa  (black).  The  MLD  shown  in  the  figure  is  the  depth  at 
which  the  temperature  becomes  0.2  °C  less  than  the  SST. 

As  in  Martin  (1985,  1986),  the  SST  is  too  high  in  summer  and,  consistent  with  this,  the  MLD 
and  ISODs  are  too  shallow.  In  Fig.  2,  the  summer  SST  is  about  2.5  °C  too  high  and  the  MLD  and 
ISODs  are  10  to  20  m  too  shallow  in  August  and  September.  The  situation  in  Fig.  2  is  a  bit  worse 
than  in  Martin  (1985,  1986)  because,  for  this  simulation,  the  net  surface  heat  flux  has  been  set 
slightly  higher  to  provide  more  net  heating  and  better  agreement  with  the  observed  upper-ocean 
heat  content  at  Papa  in  the  summer  (Martin  et  al.  2012). 
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Fig.  3  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  KC04  SDC  and  with  the 
KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 


Figure  3  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  using 
the  MYL2.5  MLM  with  the  KC04  SDC  and  with  the  KC04  LC  mixing  enhancement.  Figure  3 
shows  that  the  simulation  is  improved  relative  to  the  observations,  but  the  SST  is  still  about  1.6  °C 
too  high  in  August  and  September  and,  consistent  with  this,  the  MLD  and  ISODs  are  10  to  15  m 
too  shallow.  Hence,  the  inclusion  of  the  KC04  LC  mixing  increases  the  predicted  MLD,  but  the 
increase  does  not  seem  sufficient  to  account  for  the  observed  MLD  at  Papa  in  the  summer. 


Figure  4  shows  hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients 
Km  in  the  SML  for  the  Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement. 
Plots  of  the  ratio  of  both  unfiltered  and  temporally  filtered  Km  values  are  shown. 
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Fig.  4  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for  Papa 
simulations  with  and  without  the  KC04  LC  mixing  enhancement.  The  SDC  for  both  of  the  simulations  is  computed 
using  the  method  of  KC04.  The  left  plot  shows  the  ratio  of  unfiltered  Km  values  and  the  right  plot  shows  the  ratio 
of  temporally  filtered  Km  values. 
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The  unfiltered  plot  shows  that  the  range  of  the  ratio  of  the  Km  values  is  usually  between 
about  0.5  and  3.5,  with  occasional  lower  and  higher  values.  Values  less  than  one  indicate  that  the 
maximum  value  of  Km  for  the  simulation  with  the  KC04  LC  mixing  enhancement  is  lower  than 
that  for  the  simulation  without  the  mixing  enhancement.  This  may  sometimes  be  occurring  when 
the  vertical  shear  of  the  ocean  model’s  Eulerian  velocity  and  the  SDC  are  of  different  sign,  so  that 
the  product  of  the  two  makes  a  negative  contribution  to  the  total  shear  production  of  TKE. 

The  plot  of  the  ratio  of  the  temporally  filtered  Km  values  shows  a  range  between  about  1.0 
and  2.7.  The  ratio  is  smaller  in  the  fall  and  winter  when  free  convection,  due  mainly  to  surface 
cooling,  plays  a  larger  role  in  the  vertical  mixing,  and  is  larger  in  the  spring  and  summer  when 
forced  convection  due  to  wind  forcing  and  the  resulting  shear  production  of  TKE  plays  a  larger 
role.  The  mean  value  of  the  ratio  in  the  spring  and  summer  is  about  1.9. 

Figure  5  shows  hourly  values  of  the  magnitude  of  the  surface  current  for  the  Papa  simulations 
without  and  with  the  KC04  LC  mixing  enhancement  (note  that  these  time  series  have  been  tempo¬ 
rally  filtered).  The  mean  value  of  the  magnitude  of  the  surface  current  is  reduced  about  15%  when 
the  KC04  LC  mixing  enhancement  is  used.  This  is  somewhat  less  than  the  reduction  of  the  surface 
current  by  the  KC04  LC  mixing  enhancement  shown  in  Fig.  1.  This  is  because  of  the  frequent 
occurrence  of  inertial  oscillations  in  the  Papa  simulations,  which  tend  to  be  fairly  uniform  with 
depth  in  the  SML  and  are  not  significantly  affected  by  the  strength  of  the  vertical  mixing,  whereas 
for  the  velocity  profiles  in  Fig.  1  the  inertial  oscillations  were  suppressed  (by  turning  on  the  wind 
forcing  gradually  using  a  linear  ramp  over  one  inertial  period). 

4.2  Papa  Simulations  with  SDC  Computed  From  Buoy  Wave  Spectra 

Nine  years  of  wind  and  wave  data  from  NOAA  Buoy  46005  in  the  northeast  Pacific  were  used 
to  compute  average  wave  spectra  within  1-m/s-wide,  wind-speed  bins.  These  average  wave  spectra 
were  then  used  to  compute  hourly  SDC  profiles  at  OWS  Papa  for  1961  based  on  the  hourly  wind 
speed  at  Papa.  The  details  of  this  calculation  can  be  found  in  Appendix  C.  Buoy  46005  is  located 
at  131.001  °W,  46.100  °N,  which  is  about  1100  km  ESE  of  the  location  of  OWS  Papa.  Hence,  the 
wind  and  wave  conditions  in  the  area  of  the  buoy  should  be  roughly  similar  to  those  near  OWS 
Papa.  This  calculation  of  the  SDC  at  Papa  will  be  referred  to  as  the  Buoy  Wave  Spectra  (BWS) 
SDC. 

Figure  6  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  for 
1961  with  the  BWS  SDC,  but  without  the  KC04  LC  mixing  enhancement.  Figure  6  shows  results 
similar  to  Fig.  2,  which  used  the  KC04  SDC.  The  predicted  SST  is  about  3  °C  too  high  in  August 
and  September  and  the  predicted  MLD  at  this  time  is  about  10  to  20  m  too  shallow. 

Figure  7  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  with 
the  BWS  SDC  and  with  the  KC04  LC  mixing  enhancement.  The  predicted  SST  and  MLD  are 
slightly  improved  relative  to  Fig.6,  but  in  August  and  September,  the  SST  is  still  about  2  °C  too 
high  and  the  predicted  MLD  is  10  to  15  m  too  shallow.  The  improvement  is  slightly  less  than  for 
the  case  with  the  KC04  SDC  in  the  previous  section.  This  is  because  the  BWS  SDC  is  generally 
slightly  smaller  than  the  KC04  SDC  and  provides  a  smaller  enhancement  to  the  vertical  mixing. 

Figure  8  shows  hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients 
Km  in  the  SML  for  the  Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement. 
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Tests  of  Parameteriation  of  Langmuir- Circulation  Mixing 


Fig.  6  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  BWS  SDC,  but  without 
the  KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 
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Fig.  7  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  BWS  SDC  and  with 
the  KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 
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Fig.  8  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for 
Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement.  The  SDC  is  computed  from  the  BWS.  The 
left  plot  shows  the  ratio  of  unfiltered  Km  values  and  the  right  plot  shows  the  ratio  of  temporally  filtered  Km  values. 


Plots  of  the  ratio  of  both  the  unfiltered  and  temporally  filtered  Km  values  are  shown.  The  unfiltered 
plot  shows  the  range  of  the  ratio  of  the  Km  values  is  usually  between  about  1  and  4,  with  occasional 
lower  and  higher  values.  The  plot  of  the  ratio  of  the  temporally  filtered  Km  values  shows  a  range 
between  about  1.0  and  3.0.  As  in  Fig.  4  with  the  KC04  SDC,  the  ratio  is  larger  in  the  spring  and 
summer  when  wind-driven  forced  convection  has  a  more  dominant  role  in  the  vertical  mixing,  and 
the  mean  value  of  the  Km  ratio  in  the  spring  and  summer  is  about  2.1. 


Hence,  comparing  the  results  for  the  KC04  LC-enhanced  mixing  using  the  BWS  and  KC04 
SDCs,  they  both  increase  the  mixing  rates  in  the  SML  in  the  spring  and  summer  by  an  average 
factor  of  about  2.  However,  using  the  BWS  SDC  results  in  a  slightly  smaller  increase  in  the 
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predicted  MLD  and  a  corresponding  increase  in  the  predicted  SST  in  the  summer  because  the  SDC 
values  are  slightly  smaller. 

4.3  Papa  Simulations  with  SDC  Computed  from  SWAN  Wave  Model 

The  SWAN  wave  model  (Booij  et  al.  1999)  was  run  at  OWS  Papa  for  the  year  1961  using  the 
observed  10-m  winds  from  Papa,  and  hourly  values  of  the  SDC  were  computed  from  the  spectral 
wave  output  from  SWAN  and  used  as  input  to  NCOM.  There  are  several  advantages  to  computing 
the  SDC  from  a  wave  model  like  SWAN  over  parameterizing  the  SDC  in  terms  of  the  current  wind 
stress  as  was  done  in  the  previous  two  sections.  The  local  wave  growth  and  decay  and  the  direction 
of  the  waves  can  be  better  accounted  for  using  a  wave  model.  The  direction  of  the  waves  and  the 
SDC  can  be  somewhat  different  from  the  wind  direction,  depending  on  the  recent  history  of  the 
wind,  and  the  direction  of  the  SDC  can  vary  with  depth  since  waves  of  different  wavelengths  may 
be  propagating  in  different  directions. 

There  are  some  issues  with  regard  to  running  a  wave  model  at  a  single  location,  a  major  one 
being  how  to  parameterize  the  wave  dissipation.  Several  wave-dissipation  schemes  for  SWAN  were 
tested  by  comparison  with  wave  observations  from  NOAA  Buoy  46005,  and  a  dissipation  scheme 
was  selected  that  resulted  in  the  smallest  RMS  error  for  the  SWAN-predicted  significant  wave  height 
(this  is  discussed  in  Appendix  D).  A  related  issue  in  running  a  wave  model  at  a  single  location  is 
that  swell  propagating  from  other  areas  cannot  be  accounted  for. 

Figure  9  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  with 
the  SWAN  SDC,  but  without  the  LC  mixing  enhancement  of  KC04.  In  August  and  September, 
the  predicted  SST  is  about  3  °C  too  high  and  the  predicted  MLD  is  about  10  to  20  m  too  shallow. 

Figure  10  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  with 
the  SWAN  SDC  with  the  KC04  LC  mixing  enhancement.  The  predicted  SST  and  MLD  are  slightly 
improved.  In  August  and  September,  the  SST  is  still  about  2.5  °C  too  high  and  the  predicted  MLD 
is  10  to  20  m  too  shallow.  The  predicted  SST  in  Fig.  10  is  slightly  higher  than  for  the  corresponding 
runs  in  the  previous  two  sections  in  Fig.  3  and  7.  The  reason  for  this  is  that  the  SWAN  SDC  is, 
on  average,  smaller  than  both  the  KC04  and  BWS  SDCs. 

Figure  11  shows  hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coeffi¬ 
cients  Km  in  the  SML  for  the  Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement 
with  the  SWAN  SDC.  Plots  of  the  ratio  of  both  the  unfiltered  and  temporally  filtered  Km  values 
are  shown.  The  ratio  of  the  filtered  values  of  Km  in  Fig.  11  shows  lower  values  than  for  the  previ¬ 
ous  simulations,  i.e. ,  the  general  range  is  from  about  0.5  to  2  with  a  mean  during  the  spring  and 
summer  of  about  1.4,  whereas,  for  the  simulations  with  the  KC04  and  BWS  SDCs,  the  mean  ratio 
during  the  summer  is  about  2.  As  in  Fig.  4  and  8,  the  ratio  is  generally  larger  in  the  spring  and 
summer  than  in  the  fall  and  winter. 

As  was  discussed  in  Section  2,  the  form  of  the  SDC  shear-production  terms  in  the  TKE  and 
TLS  equations  (i.e.,  in  Eq.  (5)  and  (6))  allows  for  negative  values  if  the  vertical  gradients  of  the 
ocean-model-predicted  velocity  and  the  Stokes-drift  velocity  are  of  different  sign.  It  is  suspected 
that  this  is  part  of  the  reason  that  the  ratio  of  Km  values  for  runs  with  and  without  the  KC04 
mixing  enhancement  in  Fig.  4,  8,  and  11  sometimes  show  values  below  one.  Since  the  purpose  of 
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Fig.  9  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  SWAN  SDC,  but 
without  the  KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 
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Fig.  10  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  SWAN  SDC  and  with 
the  KC04  LC  mixing  enhancement  (red).  ISODs  are  temporally  filtered. 
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Fig.  11  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for 
Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement.  The  left  plot  shows  the  ratio  of  unfiltered 
Km  values  and  the  right  plot  shows  the  ratio  of  temporally  filtered  Km  values.  The  SDC  is  from  SWAN. 
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the  SDC  shear-production  terms  in  the  TKE  and  TLS  equations  is  to  parameterize  the  effects  of 
LC  mixing,  one  might  wonder  if  it  is  appropriate  for  these  terms  to  sometimes  reduce  the  rate  of 
mixing  in  the  SML.  To  investigate  the  extent  to  which  negative  values  of  the  SDC  shear-production 
terms  affect  the  rate  of  mixing  in  the  SML  at  Papa,  a  run  was  made  with  the  SDC  from  SWAN 
with  the  KC04  mixing  enhancement  using  the  absolute  value  of  the  SDC  shear-production  terms 
to  force  them  to  be  positive,  i.e. ,  the  SDC  shear-production  terms  in  Eq.  (5)  and  (6)  are  computed 
using 

|  |  ,  |  dv_9vs_  |  , 

dz  dz  8z  dz 


Figure  12  shows  the  Km  ratio  for  runs  with  and  without  the  KC04  mixing  enhancement  as  in 
Fig.  11  except  that,  for  the  run  with  the  KC04  mixing  enhancement,  the  SDC  shear-production 
terms  are  forced  to  be  positive.  The  result  is  that  occurrences  of  the  Km  ratio  below  one  are 
significantly  reduced,  especially  in  the  spring  and  summer.  The  ratio  of  the  temporally  filtered  Km 
values  in  Fig.  12  rarely  drops  below  one  and  the  mean  value  of  the  Km  ratio  in  the  spring  and 
summer  is  increased  slightly  from  about  1.4  to  1.6.  Since  Fig.  12  still  shows  some  values  of  the 
Km  ratio  below  one,  it  can  be  said  that  negative  values  of  the  SDC  shear-production  terms  are 
responsible  for  a  significant  fraction  of  the  occurrences  of  the  Km  ratio  being  below  one,  but  that 
there  are  some  other  causes  as  well,  which  may  just  be  due  to  differences  in  the  temporal  evolution 
of  the  mixing. 

As  noted  previously,  the  waves  and  SDCs  predicted  by  SWAN  are  frequently  not  in  the  same 
direction  as  the  wind,  and  the  direction  of  the  SDC  can  vary  with  depth  due  to  waves  of  different 
wavelengths  traveling  in  different  directions.  To  look  at  the  effect  of  this  on  the  KC04  LC  mixing 
parameterization,  a  simulation  was  conducted  at  Papa  with  the  SDC  from  SWAN  adjusted  to 
always  be  in  the  same  direction  as  the  surface  wind.  The  results  from  this  simulation  (not  shown) 
are  similar  to  the  results  shown  in  Fig.  10,  i.e.,  forcing  the  direction  of  the  SWAN  SDC  to  be 
downwind  did  not  significantly  affect  the  predicted  SST,  MLD,  ISODs,  or  the  rate  of  mixing  within 
the  SML. 


Li  et  al.  (2005)  describe  a  regime  diagram  for  classifying  turbulent  mixing  in  the  upper  ocean. 
Such  a  diagram  is  presented  here  in  Fig.  13  for  OWS  Papa  for  1961  using  the  SDC  from  SWAN. 
The  figure  shows  a  scatter  plot  of  hourly  values  of  the  Hoenikker  number  versus  the  turbulent 
Langmuir  number  for  Papa  for  1961.  The  Hoenikker  number  H0  is  given  by 


4:B0ds 

Fs(0)ti*  ’ 


(8) 


where  Bq  is  the  total  surface  buoyancy  flux,  ds  is  the  e-folding  depth  of  the  SDC,  14(0)  is  the  surface 
value  of  the  SDC,  and  w*  is  the  water-side  surface  friction  velocity.  The  turbulent  Langmuir  number 
Leif  is  given  by 


Lat  = 


(9) 


The  diagram  is  divided  into  three  regimes  according  to  the  type  of  turbulent  mixing:  Langmuir, 
convection,  or  shear,  that  dominates  within  that  region.  The  boundaries  of  the  three  regimes 
were  determined  by  Li  et  al.  (2005)  based  on  LES  simulations.  The  regimes  are  characterized  by 
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Fig.  12  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for 
Papa  simulations  with  and  without  the  KC04  LC  mixing  enhancement.  The  left  plot  shows  the  ratio  of  unfiltered 
Km  values  and  the  right  plot  shows  the  ratio  of  temporally  filtered  Km  values.  For  the  simulation  with  the  KC04 
mixing  enhancement,  the  SDC  shear-production  terms  are  forced  to  be  positive. 
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Hoen i kker  Number  vs  Turbulent  Langmuir  Number 
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LoglO  (Turbulent  Langmu i n  Numben) 

Fig.  13  —  Turbulent  regime  diagram  for  OWS  Papa  for  1961  with  SDC  from  SWAN.  The  points  plotted  are  hourly 
values  of  the  loglO  of  the  Hoenikker  number  (with  an  offset  of  0.01)  versus  the  loglO  of  the  turbulent  Langmuir 
number. 
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the  relative  magnitudes  of  the  components  of  the  turbulent  mixing,  with  the  vertical  component 
dominating  in  convective  mixing,  the  downwind  component  dominating  in  shear  mixing,  and  the 
cross-wind  and  vertical  components  dominating  the  downwind  component  in  Langmuir  mixing. 

Figure  13  indicates  that  the  upper-ocean  mixing  at  Papa  for  1961  is  generally  in  the  Langmuir 
and  convection  regimes.  Only  a  few  points  occur  in  the  shear  regime,  and  these  occur  close  to  the 
boundary  between  the  shear  and  Langmuir  regimes.  This  is  consistent  with  the  finding  of  Li  et  al. 
(2005)  that  wind  mixing  in  the  upper-ocean  typically  occurs  in  the  Langmuir  regime. 

Regime  diagrams  for  Papa  for  1961  were  also  calculated  (not  shown)  with  the  SDC  computed 
using  (a)  the  formulas  from  KC04  and  (b)  the  mean  wave  spectra  at  different  wind  speeds  from 
Buoy  46005  (as  discussed  in  the  previous  two  sections  and  in  Appendices  B  and  C,  respectively). 
These  plots  are  similar  to  Fig.  13,  though  the  range  of  La./  is  smaller  since  the  waves  and  their 
associated  SDC  are  computed  based  on  the  current  wind  speed.  For  both  of  these  cases,  Leif  has  a 
value  near  0.3,  which  is  close  to  the  mean  value  for  the  SWAN  SDC  in  Fig.  13,  and  is  also  close  to 
the  value  for  typical  ocean  conditions  reported  by  Li  et  al.  (2005). 

5.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  FOR  HURRICANE  IVAN 

Simulations  were  run  with  NCOM  in  the  Gulf  of  Mexico  (GoM)  with  wave  forcing  from  SWAN 
to  look  at  the  effect  of  the  KC04  LC  mixing  enhancement  on  a  simulation  of  the  response  of  the 
ocean  to  Hurricane  Ivan  in  September  2004.  Atmospheric  forcing  was  provided  from  a  simulation 
of  Hurricane  Ivan  with  the  COAMPS  atmospheric  model.  COAMPS,  which  stands  for  Coupled 
Ocean/ Atmosphere  Mesoscale  Prediction  System  (Hodur  1997),  is  currently  a  coupled  atmosphere- 
ocean-wave  modeling  system.  However,  for  the  simulations  conducted  here,  the  atmosphere,  ocean, 
and  wave  models  were  run  separately.  Within  the  rest  of  this  report,  COAMPS  will  be  used  to 
refer  to  just  the  COAMPS  atmospheric  model. 

NCOM  was  run  with  a  horizontal  grid  resolution  of  about  4  km.  NCOM’s  vertical  grid  con¬ 
sisted  of  40  layers  with  a  surface  layer  thickness  of  1  m  and  a  smooth  stretching  to  a  maximum 
depth  of  5500  m.  Initial  and  boundary  conditions  for  NCOM  were  provided  by  the  Global  NCOM 
model  that  is  run  operationally  at  the  Naval  Oceanographic  Office  at  Stennis  Space  Center,  Mis¬ 
sissippi  (Barron  et  al.  2004).  The  COAMPS  atmospheric  model  was  run  as  a  triply- nested  system 
with  grid  resolutions  of  81,  27,  and  9  km.  Initial  and  boundary  conditions  for  COAMPS  were 
provided  by  the  Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS),  which  is 
run  at  the  Navy’s  Fleet  Numerical  Meteorology  and  Oceanography  Center  (FNMOC)  in  Monterey, 
California  (Rosmond  et  al.  2002).  Hourly  surface  fields  from  COAMPS’s  innermost  grid  were  used 
to  force  SWAN  and  NCOM.  These  fields  include  the  surface  winds  for  SWAN,  and  the  surface  wind 
stress,  surface  pressure,  solar  radiation,  net  longwave  radiation,  latent,  and  sensible  heat  fluxes, 
evaporation,  and  precipitation  for  NCOM.  SWAN  was  run  on  the  same  horizontal  grid  as  NCOM. 

Wave  fields  for  forcing  NCOM  were  output  hourly  from  SWAN  and  include  the  surface  wave 
radiation  stress,  the  SDC  profiles,  and  the  characteristic  amplitude,  frequency,  and  direction  of 
the  orbital  wave  motion  near  the  ocean  bottom.  The  wave  radiation  stress  forcing  in  NCOM  is 
implemented  similar  to  the  surface  wind  stress  as  a  surface  stress.  In  deep  water,  the  wave  radiation 
stress  forcing  of  the  ocean  is  relatively  small,  i.e.,  generally  less  than  10%  of  the  magnitude  of  the 
wind  stress.  The  SDC,  as  noted  in  Section  2,  is  used  in  NCOM’s  Coriolis  and  advection  terms,  in 
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the  continuity  equation,  and  also  in  the  KC04  mixing  enhancement.  The  wave  orbital  motions  near 
the  bottom  act  to  enhance  the  ocean  model’s  bottom  drag  in  areas  where  the  bottom  is  sufficiently 
shallow  to  feel  the  orbital  wave  motions  (Grant  and  Masden  1979;  Soulsby  1995).  However,  the 
focus  here  is  only  on  the  effect  of  the  KC04  mixing  enhancement;  hence,  the  two  NCOM  simulations 
that  are  compared  both  have  full  wave  forcing  from  SWAN,  and  the  only  difference  between  them 
is  whether  the  KC04  mixing  enhancement  is  used  or  not. 

The  ocean  model  simulation  was  started  at  00Z  on  September  12,  2004.  Figure  14  shows  the 
SST  and  sea-surface  velocity  (SSV)  vectors  for  the  simulation  with  the  KC04  mixing  enhancement 
for  09Z  September  15.  The  surface  wind  vectors  from  the  COAMPS  atmospheric  model  are  also 
shown.  Note  that  the  surface  current  vectors  in  Fig.  14  include  both  the  NCOM  surface  current  and 
the  SDC  from  SWAN.  Hurricane  Ivan  is  near  the  center  of  the  GoM  at  this  time.  The  hurricane 
winds  generate  strong  cooling  along  the  hurricane  track,  with  the  strongest  cooling  usually  occurring 
just  to  the  right  of  the  hurricane  track  (Price  1981). 

Figure  15  shows  the  difference  in  the  SST  between  the  simulation  with  the  KC04  mixing 
enhancement  and  without  at  09Z  September  15,  2004.  The  simulation  with  the  KC04  mixing 
enhancement  generally  shows  slightly  more  cooling  of  the  SST  as  indicated  by  the  areas  of  negative 
temperature  difference,  but  the  difference  is  generally  small  (less  than  0.5  °C).  The  slightly  lower 
SST  for  the  simulation  with  the  KC04  mixing  enhancement  indicates  slightly  deeper  mixing,  which 
is  consistent  with  the  previous  results  in  Sections  3  and  4. 

Figure  16  shows  the  difference  in  the  SSV  between  the  simulation  with  the  KC04  mixing 
enhancement  and  without  at  09Z  September  15,  2004.  The  color  contours  show  the  difference  in 
the  magnitude  of  the  surface  current  and  the  black  vectors  show  the  vector  velocity  difference.  The 
fact  that  the  direction  of  the  velocity  difference  vectors  is  opposite  to  the  wind  direction  indicates 
that  the  stronger  mixing  with  the  KC04  mixing  enhancement  reduces  the  wind-driven  surface 
currents.  Near  the  hurricane,  the  surface  currents  are  reduced  by  20  cm/s  or  more.  Some  areas 
show  the  surface  current  magnitude  with  the  KC04  mixing  enhancement  to  be  larger  (as  indicated 
by  the  color  contours  in  Fig.  16)  and  this  is  because  the  direction  of  the  prevailing  currents  in  these 
areas  is  nearly  opposite  to  the  wind  direction,  and  so  a  reduction  in  the  downwind,  wind-driven 
current  increases  the  total  current,  e.g.,  on  the  west  side  of  the  hurricane  there  is  an  anti-cylonic 
Loop  Current  eddy  (see  Fig.  14). 

Figure  17  shows  the  maximum  value  of  the  vertical  mixing  coefficient  Km  in  the  SML  for  the 
simulation  with  the  KC04  mixing  enhancement  at  09Z  September  15,  2004.  The  Km  values  in 
Fig.  17  were  spatially  filtered  with  one  pass  of  a  3x3  box  filter.  The  strongest  mixing  near  the 
hurricane  is  occurring  on  the  north  and  east  sides  of  the  eye.  A  circular  area  of  weak  mixing  occurs 
under  the  eye  of  the  hurricane  due  to  the  weak  winds  there. 

The  local  time  for  Fig.  17  is  about  three  o’clock  in  the  morning;  hence,  there  is  relatively  strong 
night-time  convective  mixing  occurring  in  some  areas,  e.g.,  in  the  western  GoM  and  the  northwest 
Caribbean.  On  the  afternoon  of  September  15  (local  time),  most  of  this  mixing  is  significantly 
reduced  due  to  day-time  solar  heating  (not  shown).  However,  the  mixing  near  the  hurricane  is 
primarily  wind-driven  and  remains  strong  through  the  day  and  night.  Also,  the  stronger  mixing 
south  of  and  just  north  of  Cuba  is  partly  due  to  this  area  being  on  the  warm  side  of  the  Loop 
Current  where  the  SML  is  generally  deeper  than  on  the  cold  side. 
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Fig.  14  —  SST  (color  contours),  SSV  (black  vectors),  and  surface  winds  (white  vectors)  at  09Z  September  15,  2004 
for  simulation  of  Hurricane  Ivan  with  the  KC04  mixing  enhancement 
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Fig.  15  —  Difference  in  SST  for  Hurricane  Ivan  simulations  with  and  without  the  KC04  mixing  enhancement  at  09Z 

September  15,  2004. 
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Fig.  16  —  Difference  in  SSV  for  Hurricane  Ivan  simulations  with  and  without  the  KC04  mixing  enhancement  at  09Z 
September  15,  2004.  The  color  contours  show  the  difference  in  the  magnitude  of  the  surface  current  and  the  black 
vectors  show  the  vector  velocity  difference.  The  white  vectors  show  the  surface  winds. 
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Fig.  17  —  Maximum  value  of  Km  in  SML  for  Hurricane  Ivan  simulations  with  KC04  mixing  enhancement  at  09Z 

September  15,  2004.  The  white  vectors  show  the  surface  winds. 
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Fig.  18  —  Ratio  of  maximum  values  of  Km  in  SML  for  Hurricane  Ivan  simulations  with  and  without  the  KC04  mixing 
enhancement  at  09Z  September  15,  2004.  The  white  vectors  show  the  surface  winds. 


Figure  18  shows  the  ratio  of  the  spatially-filtered  values  of  the  maximum  value  of  Km  in  the 
SML  for  the  simulations  with  and  without  the  the  KC04  mixing  enhancement  at  09Z  September 
15,  2004.  Near  the  hurricane,  the  ratios  range  from  1  to  4  and  are  frequently  in  the  range  of  2  to 
3.  This  indicates  that  the  KC04  scheme  is  significantly  enhancing  the  mixing  rates  in  the  areas  of 
strong  hurricane  winds.  The  areas  where  the  ratio  is  below  one  are  generally  occurring  in  areas  of 
night-time  convection. 

6.  SUMMARY 

Mellor-Yamada-type  turbulence  models  have  long  been  criticised  for  not  mixing  the  ocean’s 
surface  layer  strongly  or  deeply  enough  in  the  open  ocean.  Recent  LES  simulations  of  Langmuir 
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circulation  in  the  upper  ocean,  e.g.,  MW97,  found  that  the  occurrence  of  LCs  significantly  increases 
the  net  rate  of  mixing  within  the  SML  and  slightly  increases  the  MLD.  Mellor-Yamada-type  tur¬ 
bulence  models  do  not  currently  account  for  the  effect  of  LCs;  hence,  based  on  the  LES  results  of 
MY97,  KC04  parameterized  and  tested  the  effects  of  LC  mixing  in  the  frequently  used  MYL2.5 
turbulence  model. 

KC04  parameterized  the  effect  of  LC  mixing  in  the  MYL2.5  turbulence  model  by  adding  ad¬ 
ditional  shear-production  terms  to  the  TKE  and  TLS  equations  of  the  turbulence  model.  These 
additional  terms  consist  of  the  product  of  the  vertical  turbulent  momentum  flux  and  the  vertical 
shear  of  the  SDC  from  the  surface  wave  field. 

The  KC04  LC  mixing  parameterization  was  implemented  in  the  MYL2.5  turbulence  model  in 
NCOM  and  tested  by  conducting  simulations  of  the  simple,  light-wind,  mixing  case  used  by  both 
MW97  and  KC04,  the  upper-ocean  thermal  structure  at  OWS  Papa  during  1961,  which  was  also 
used  by  KC04,  and  Hurrican  Ivan  in  the  Gulf  of  Mexico  in  2004. 

Our  initial  implementation  of  the  KC04  LC  mixing  parameterization  did  not  perform  as  well 
as  reported  by  KC04,  and  this  turned  out  to  be  due  to  an  error  in  the  value  of  the  constant  used 
to  scale  the  SDC  shear-production  term  in  the  TLS  equation.  The  value  of  the  constant  had 
been  incorrectly  reported  in  the  KC04  paper  as  4.0,  and  we  were  advised  by  Dr.  Lakshmi  Kantha 
(personal  communication)  that  it  should  have  a  value  of  7.2. 

Our  simulation  of  the  simple  test  case  used  by  MW97  and  KC04,  which  consists  of  a  constant 
surface  wind  stress  of  0.037  Pa  (which  corresponds  to  a  wind  speed  of  about  5  m/s)  and  an  initial 
MLD  of  33  m,  gave  results  similar  to  those  reported  by  both  MW97  and  KC04.  The  addition  of 
the  KC04  LC  mixing  parameterization  increased  the  rate  of  mixing  by  about  a  factor  of  two  and 
increased  the  MLD  by  about  3  m.  MW97  actually  showed  an  increase  in  the  rate  of  mixing  of 
about  a  factor  of  three,  but  the  maximum  mixing  rate  in  the  SML  for  all  three  simulations  with 
the  LC  mixing  was  similar,  i.e.,  about  400  m2/s. 

For  the  Papa  simulations  for  1961,  no  wave  observations  are  available.  Hence,  we  used  three 
methods  to  estimate  the  SDC:  (1)  the  parameterization  of  the  SDC  in  terms  of  the  surface  friction 
velocity  u*  used  by  KC04  in  their  tests  at  Papa,  (2)  the  calculation  of  the  mean  wave-energy 
spectrum  at  different  wind  speeds  using  data  from  NOAA  Buoy  46005  in  the  northeast  Pacific, 
and  (3)  the  waves  predicted  by  the  SWAN  wave  model,  which  was  run  using  the  observed  winds  at 
Papa  for  1961. 

For  all  three  calculations  of  the  SDC  at  Papa,  the  use  of  the  KC04  LC  mixing  parameterization 
significantly  increased  the  rate  of  mixing  within  the  SML  and  increased  the  MLD.  However,  the 
methods  (1)  to  (3)  of  computing  the  SDC  yielded  sucessively  smaller  values  of  the  SDC,  at  least 
during  the  heating  season  in  the  spring  and  summer,  which  resulted  in  slightly  reduced  mixing 
rates  and  reduced  deepening  of  the  SML.  The  mean  rate  of  mixing  within  the  SML  in  the  spring 
and  summer  was  increased  by  a  factor  of  about  2  for  methods  (1)  and  (2)  and  by  a  factor  of  about 
1.4  for  method  (3).  Note  that  for  method  (1),  KC04  indicated  that  the  scaling  parameter  they  used 
to  compute  the  SDC  from  u *  was  on  the  high  side;  hence,  the  fact  that  the  SDCs  computed  in  this 
way  for  OWS  Papa  for  1961  are  significantly  higher  than  the  SWAN  SDCs  suggests  that  this  is  the 
case. 
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The  use  of  the  KC04  LC  mixing  parameterization,  besides  increasing  the  rate  of  mixing  within 
the  SML,  also  increases  the  MLD.  However,  for  the  basic  MYL2.5  turbulence  model,  the  increase 
in  the  MLD  provided  by  the  KC04  LC  mixing  parameterization  is  not  deep  enough  to  match  the 
observed  MLD  at  Papa  for  1961.  Hence,  it  may  be  that  an  additional  mechanism  is  needed  to 
provide  the  needed  additional  deepening,  e.g.,  the  parameterization  of  intermittent  mixing  below 
the  SML  (e.g.,  by  internal  waves)  of  Large  et  al.  (1994). 

The  additional  shear-production  terms  in  the  TKE  and  TLS  equations  implemented  by  KC04, 
which  consist  of  the  product  of  the  vertical  momentum  flux  and  the  vertical  shear  of  the  SDC, 
can  make  a  negative  contribution  to  the  total  shear  production  of  TKE  if  the  vertical  shear  of  the 
ocean  model  velocity  and  the  vertical  shear  of  the  SDC  are  of  opposite  sign.  The  fact  that  mixing 
rates  with  the  KC04  mixing  enhancement  are  sometimes  less  than  without  it  suggests  that  this 
sometimes  occurs.  A  simulation  was  conducted  at  Papa  with  the  SWAN  SDC  with  KC04’s  SDC- 
shear-production  terms  in  the  TKE  and  TLS  equations  forced  to  make  a  positive  contribution  to 
the  total  shear  production  of  TKE  by  using  their  absolute  value.  The  result  was  that  the  occurrence 
of  reduced  mixing  rates  when  using  the  KC04  LC  mixing  parameterization  was  reduced,  and  the 
ratio  of  the  mean  value  of  the  vertical  mixing  coefficients  during  the  spring  and  summer  for  the 
cases  with  and  without  the  LC-mixing  parameterization  was  increased  slightly  from  1.4  to  1.6. 

In  computing  the  SDC  at  Papa  using  methods  (1)  and  (2),  the  SDC  was  taken  to  be  in  the 
same  direction  as  the  wind.  However,  the  waves  predicted  by  SWAN  in  method  (3)  (and,  of  course, 
real  waves)  often  propagate  in  a  direction  different  from  that  of  the  wind,  and  waves  of  different 
wavelengths  may  propagate  in  different  directions,  which  results  in  a  change  in  direction  of  the 
SDC  with  depth.  The  effect  of  this  on  the  contribution  to  the  vertical  mixing  from  the  KC04 
LC  mixing  parameterization  was  investigated  by  conducting  a  simulation  at  Papa  with  the  SWAN 
SDC  modified  to  always  be  alligned  with  the  wind.  The  effect  of  this  change  on  the  simulation  was 
small,  i.e.,  the  SST,  MLD,  and  mixing  rates  within  the  SML  were  not  much  affected. 

For  the  simulations  of  Hurrican  Ivan  in  the  Gulf  of  Mexico,  the  use  of  the  KC04  LC  mixing 
parameterization  generally  results  in  significantly  increased  mixing  rates  in  the  SML  near  the 
hurricane,  i.e.,  the  mixing  rates  are  frequently  increased  by  a  factor  of  2  to  3.  The  increased  mixing 
rates  result  in  a  significant  reduction  in  the  vertical  shear  of  the  wind-driven  currents  in  the  SML, 
which  results  in  a  reduction  of  the  surface  value  of  the  wind-driven  current  near  the  hurricane  of 
20  cm/s  or  more.  The  depth  of  mixing  is  also  slightly  increased,  which  results  in  increased  cooling 
of  the  SST,  though  this  additional  SST  cooling  is  small,  i.e.,  generally  less  than  0.5  °C. 
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CALCULATION  OF  SDC  FOR  A  MONOCROMATIC  WAVE 


The  Stokes  drift  current  (SDC)  Vs  can  be  computed  for  a  monocromatic  ( i.e. ,  single  frequency) 


wave  using 


V.(z)  =  (aO%COSh|Mf  +  211. 

V  ’  K  p  2sinh 2{kH) 


where  a  is  the  wave  amplitude,  k  =  2-7r/L  is  the  wavenumber,  L  is  the  wavelength,  cp  = 
(t&nh(kH)g/k)z  is  the  phase  speed,  g  is  the  acceleration  of  gravity,  H  is  the  total  water  depth, 
and  z  is  the  distance  to  the  surface. 


In  deep  water,  this  simplifies  to 

Vs(z)  =  (ak)2cpex^p(2kz),  (A2) 

where  the  phase  speed  is  cp  =  (g/k)  t, 

For  a  wave  field  consisting  of  waves  of  different  amplitudes,  frequencies,  and  directions,  the 
contributions  to  the  SDC  from  the  different  waves  can  be  summed  (in  a  vector  sense). 

The  wave  amplitude  a  is  related  to  the  wave  energy  E  as 

E  =  in2.  (A3) 

Note  that  the  wave  energy  given  by  Eq.  (A3)  has  units  of  length  squared,  which  is  a  common 
practice  when  discussing  surface  waves.  Multiplication  by  the  water  density  and  the  acceleration 
of  gravity  gives  the  proper  units  of  energy  per  m2  of  surface  area. 
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CALCULATION  OF  SDC  FROM  WIND  STRESS  BY  KC04 

Kantha  and  Clayson  (2004)  estimated  the  surface  value  of  the  SDC  Vs(z  =  0)  in  deep  water 
using 


Vs(z  =  0)  =  11.8t»*,  (Bl) 

where  u*  =  {r / pw)  2  is  the  water-side  surface  friction  velocity,  r  is  the  surface  wind  stress,  and  pw 
is  the  water  density.  The  wavenumber  k  associated  with  the  peak  of  the  wave  energy  spectrum  is 
parameterized  as 


k  =  4.05  x  10“  V^. 


Hence,  the  SDC  profile  is  given  by 


Vs(z)  =  11.8u*exp(2kz), 


with  k  given  by  Eq.  B2. 


(B2) 


(B3) 
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CALCULATION  OF  SDC  USING  WAVE  SPECTRA  FROM  NOAA  BUOY 


Some  of  the  NOAA  buoys  provide  observations  of  both  the  surface  wind  velocity  and  the  wave 
energy  spectra.  These  data  allow  the  calculation  of  the  average  wave  spectra  at  different  wind 
speeds.  The  average  wave  spectra  at  the  different  wind  speeds  can  then  be  used  to  estimate  the 
SDC  for  a  given  wind  speed. 

The  spectral  wind  and  wave  data  that  were  used  were  from  NOAA  Buoy  46005,  which  is  located 
in  the  northeast  Pacific  at  131.001  °W,  46.100  °N,  which  is  about  1100  km  east-southeast  of  the 
location  of  OWS  Papa.  Since  the  location  of  the  buoy  is  somewhat  similar  to  that  of  OWS  Papa, 
the  wind  and  wave  conditions  should  be  roughly  similar. 

Data  from  Buoy  46005  for  the  years  1996  through  2004  were  used.  For  these  years,  the  spectral 
wave  data  from  the  buoy  consist  of  38  values  of  the  wave  energy  at  frequencies  from  0.03  to  0.4  Hz 
at  0.01  Hz  intervals.  Note  that  no  directional  information  is  provided  for  the  wave  data  from  Buoy 
46005.  Averages  of  the  spectral  data  were  computed  within  1-m/s-wide  wind-speed  bins  from  1 
to  20  m/s.  No  winds  above  20  m/s  were  recorded  at  Buoy  46005  during  the  9-year  period.  The 
wind  speed  used  to  characterize  the  wave  spectrum  in  this  calculation  is  the  wind  speed  at  the  time 
of  the  wave  observation.  An  alternative  would  be  to  use  some  time-weighted  wind  speed  over  the 
previous  few  hours,  but  this  was  not  done. 

Figure  Cl  shows  the  averaged  wave  spectra  from  NOAA  Buoy  46005  for  wind  speeds  of  1,  5, 
10,  15,  and  20  m/s.  The  wave  energy  increases  with  increasing  wind  speed.  There  are,  however,  a 
couple  of  notable  aspects  to  the  wave  spectra  in  Fig.  Cl:  (1)  the  location  of  the  peak  of  the  wave 
energy  spectra  shows  little  change  with  wind  speed,  and  (2)  there  is  a  notable  amount  of  wave 
energy  present  at  low  wind  speeds  of  1  to  5  m/s. 

There  are  several  limitations  to  estimating  the  SDC  from  a  wave  spectra  based  on  the  current 
local  wind  velocity.  The  wind  duration  and  fetch  and  the  propagation  of  waves  from  other  areas 
cannot  be  accounted  for.  Additionally,  variation  in  the  direction  of  the  waves  making  up  the  energy 
spectrum  cannot  be  accounted  for;  hence,  the  direction  of  the  waves  and  the  SDC  is  assumed  to 
be  downwind. 

The  SDC  profiles  estimated  from  the  wave  spectra  at  Buoy  46005  for  different  wind  speeds  were 
compared  with  the  SDC  profiles  computed  from  the  observed  wave  spectra  using  data  for  the  year 
2000,  and  the  mean  and  mis  differences  versus  depth  are  shown  in  Fig.  C2.  The  SDC  estimated 
by  this  method  will  be  referred  to  as  the  Buoy  Wave  Spectra  (BWS)  SDC.  Mean  and  RMS  errors 
with  respect  to  the  SDC  computed  from  the  buoy  wave  observations  for  2000  are  also  shown  for 
the  SDC  computed  using  the  KC04  method  discussed  in  Appendix  C2. 
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Wave  Energy 


Fig.  Cl  —  Averaged  wave  spectra  from  NOAA  Buoy  46005  for  wind  speeds  of  1,  5,  10,  15,  and  20  m/s.  The  wave 

energy  increases  with  wind  speed. 
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Depth  (m) 


Fig.  C2  —  Mean  and  RMS  error  for  the  KC04-computed  (solid  line)  and  BWS-computed  (dashed  line)  SDC  at  buoy 
46005  for  the  year  2000.  Errors  are  with  respect  to  the  SDC  computed  from  the  observed  wave  spectra  at  buoy  46005 
for  the  year  2000. 


Martin  et  al. 


Figure  C2  shows  that  both  the  mean  and  RMS  errors  are  noticeably  smaller  for  the  BWS  SDC 
than  for  the  SDC  estimated  using  the  method  of  KC04.  Plots  for  the  same  errors  computed  using 
data  from  buoy  46005  for  other  years  between  1996  and  2004  (not  shown)  look  similar  to  Fig.C2. 
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CALCULATION  OF  SDC  WITH  SWAN  WAVE  MODEL 


This  section  discusses  the  use  of  the  SWAN  wave  model  (Booij  et  al.  1999)  to  estimate  the 
wave  spectra  and  SDC  at  OWS  Papa.  For  the  application  of  SWAN  at  a  single  location,  the 
geographic  advection  terms  and  a  feature  known  as  “quadrant  sweeping”  are  disabled.  The  use  of 
a  wave  model  at  a  single  location  allows  accounting  for  the  recent  time  history  of  the  local  winds, 
including  changes  in  both  magnitude  and  direction.  However,  it  does  not  allow  a  proper  accounting 
for  swell  from  other  areas,  where  the  winds  are  likely  to  be  different. 

Rather  than  apply  the  model  with  the  Papa  winds  blindly,  we  first  perform  several  experiments 
for  a  situation  where  simultaneous  observations  of  wind  speed  and  wave  spectra  are  available,  i.e., 
with  the  data  from  NOAA  buoy  46005  during  the  year  2000.  As  noted  previously,  buoy  46005  is 
located  about  1100  km  east-southeast  of  the  location  of  OWS  Papa,  so  the  winds  should  be  roughly 
similar  in  a  statistical  sense.  With  the  buoy  data,  we  can  check  the  accuracy  of  our  method,  and 
evaluate  simulations  with  several,  alternative,  wave-model  source-term  settings.  To  simplify  the 
analysis  of  the  results,  we  look  at  the  significant  wave  height  (SWH)  only,  rather  than  the  wave 
spectra.  (An  alternative  would  be  to  look  at  an  omni-directional  calculation  of  the  Stokes  drift  at 
a  selected  depth.) 

The  physics  packages  available  with  the  version  of  SWAN  that  we  have  that  are  compared  here 
are:  (1)  the  Komen  et  al.  (1984)  physics  (denoted  here  as  KHH,  this  is  the  SWAN  default)  with 
the  Wu  (1980)  drag  formulation,  (2)  the  KHH  physics  with  modification  according  to  Rogers  et  al. 
(2003)  with  the  Hwang  (2011)  drag  formulation,  (3)  the  Rogers  et  al.  (2012)  physics  (denoted  as 
RBW)  with  the  Hwang  (2011)  drag  formulation  without  swell  dissipation,  and  (4)  variants  of  (3) 
with  the  swell  dissipation  strength  controlled  by  the  swell  dissipation  coefficient  fe.  For  further 
description  of  these  physics  packages,  the  reader  is  referred  to  Rogers  et  al.  (2012). 

The  duration  of  the  test  simulations  is  one  year.  The  mean  and  rms  errors  for  the  SWH  for  the 
tests  that  were  run  are  listed  in  Table  Dl.  Errors  for  the  previous  KC04  and  BWS  parameterizations 
of  the  wave  field  in  terms  of  the  local  winds  are  also  provided  in  Table  Dl  for  comparison. 

The  minimum  error  in  Table  Dl  is  for  the  SWAN  simulation  with  the  Rogers  et  al.  (2012) 
physics  with  the  Hwang  (2011)  drag  formulation  with  a  swell  dissipation  coefficient  of  0.003.  The 
SWH  RMS  error  for  this  simulation  is  0.99  m.  Figure  Dl  shows  a  scatter  plot  of  the  SWH  predicted 
by  SWAN  for  this  case  versus  the  observed  SWH  at  buoy  46005  during  the  year  2000. 

Though  the  SWH  RMS  error  of  0.99  m  is  the  lowest  that  was  obtained  in  the  tests  conducted 
with  the  buoy  data,  this  is  significantly  greater  than  the  SWH  RMS  error  of  0.4  to  0.6  m  typically 
expected  for  global  wave  hindcasts  (e.g.,  Bidlot  et  al.  2002).  This  larger  error  for  the  wave  simulation 
at  buoy  46005  can  be  blamed  on  the  use  of  a  single-point  wave  model,  which  has  these  limitations: 
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Table  D1  —  Significant  Wave  Height  Errors  Using  Data 
From  Buoy  46005 


model 

physics 

bias  error 
(m) 

rms  error 

(m) 

SWAN 

KHH/Wu,  n  =  1 

-1.09 

1.43 

SWAN 

KHH/Hwang,  n  =  2 

2.16 

2.37 

SWAN 

RBW,  fe  =  0.000 

9.58 

9.81 

SWAN 

RBW,  fe  =  0.002 

0.51 

1.20 

SWAN 

RBW,  fe  =  0.003 

-0.03 

0.99 

SWAN 

RBW,  fe  =  0.004 

-0.34 

1.00 

SWAN 

RBW,  fe  =  0.005 

-0.69 

1.16 

KC04 

-1.20 

1.74 

BWS 

0.28 

1.13 

(1)  it  assumes  uniform  winds  within  an  infinite  ocean  domain,  i.e.,  the  fetch  is  effectively  infinite, 

(2)  it  cannot  use  advection  to  remove  swell  energy,  the  energy  can  be  removed  only  via  dissipation, 
and  (3)  it  cannot  use  advection  to  bring  in  swell  energy,  energy  is  created  only  via  the  wind  input 
physics, 

Intuition  tells  us  that  a  single-point  wave  model  will  tend  to  do  best  during  periods  when 
the  sea  state  is  dominated  by  local  winds.  Further,  it  is  intuitive  that  models  with  weak  swell 
dissipation,  e.g.,  the  SWAN  simulations  with  the  KHH/Hwang  physics  with  n  =  2  and  with  the 
RBW  physics  with  fe  =  0,  will  tend  to  perform  poorly. 

Based  on  these  results,  the  SWAN  wave  model  with  the  RBW  physics  with  fe  =  0.003  is 
applied  to  the  Papa  wind  data  set  to  estimate  the  wave  spectra  and,  from  the  wave  spectra,  the 
Stokes  drift. 
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Fig.  D1  —  Scatter  plot  of  SWH  predicted  by  SWAN  versus  observed  SWH  at  Buoy  46005  for  year  2000  for  SWAN 
simulation  with  minimum  error,  i.e.,  with  RBW  physics  with  fe  =  0.003. 
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